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Monte Carlo simulations of a spherical macroion, surrounded by a size-asymmetric electrolyte 
in the primitive model, were performed. We considered 1:1 and 2:2 salts with a size ratio of 2 
(i.e., with coions twice the size of counterions) , for several surface charge densities of the macro- 
sphere. The radial distribution functions, electrostatic potential at the Helmholtz surfaces, and 
integrated charge are reported. We compare these simulational data with original results obtained 
from the Ornstein-Zernike integral equation, supplemented by the hypernetted chain/hypernetted 
chain (HNC/HNC) and hypernetted chain/mean spherical approximation (HNC/MSA) closures, and 
with the corresponding calculations using the modified Gouy-Chapman and unequal-radius modi- 
fied Gouy-Chapman theories. The HNC/HNC and HNC/MSA integral equations formalisms show 
good concordance with Monte Carlo "experiments", whereas the notable limitations of point- ion 
approaches are evidenced. Most importantly, the simulations confirm our previous theoretical pre- 
dictions of the non-dominance of the counterions in the size-asymmetric spherical electrical double 
layer [J. Chem. Phys. 123, 034703 (2005)], the appearance of anomalous curvatures at the outer 
Helmholtz plane and the enhancement of charge reversal and screening at high colloidal surface 
charge densities due to the ionic size asymmetry. 

PACS numbers: 61.20.Gy,61.20. Ja,61.20.Ne,61.20.Qg. 



I. INTRODUCTION 



The study of charged colloidal solutions is very relevant 
for both basic research and technology due to the ubiq- 
uitous nature of these systems [E S S B @ • Accord- 
ingly, the attainment of a successful theoretical descrip- 
tion of such state of matter should represent a keystone 
for later developments in colloid science. For many years, 
the scientific community has investigated the structural 
characteristics of these materials, trying to understand 
the role of the electrostatic and entropic correlations in 
their observable properties. In particular, the interest in 
charged suspensions has prompted the burgeoning of un- 
precedented experimental techniques and of numeric and 
statistical mechanics approaches of increasing complex- 
ity. On the theoretical side, and despite of the notorious 
progress in the speed of machine calculations, at present, 
it is not yet possible to mimic a real dispersion without 
making several and important simplifications in order to 
establish a tractable problem. Thus, for example, one 
of the most elemental idealizations of a diluted charged 
colloidal suspension is the combination of the cell and 
primitive models. Within this scheme, the average dis- 
tance between non-concentrated macroions bathed by an 
electrolyte is very large, and therefore it is expected that 
the thermodynamics of the system will depend mainly 
on the ionic structure, or electrical double layer (EDL), 
around a single macroparticle enclosed in an electroneu- 
tral cell. Complementarily, the so-called primitive model 
(PM), in which the ions are treated as hard spheres with 
punctual charges embedded in their centers and the sol- 
vent is considered a continuous medium, stands as the 
most thriving representation of a multi-component elec- 



trolyte. A particular case of the PM is the restricted 
primitive model (RPM) , where all the ionic species are of 
equal size. This condition drastically facilitates the the- 
oretical analysis and, as a consequence, a great amount 
of work has been performed in the RPM for the planar 
0SJS[S[pJllL cylindrical [II QUE El and spher- 
ical [Hill HaEa [ID, [H SH double layers. In strong 
contrast, there are few articles in which the effects of ionic 
size asymmetry have been studied systematically, and 
the se p ublications focus chiefly on the planar instance 

Sllal^l^l2El2il^|3il^|3lS^I^. 



Certainly, the widespread use of the RPM to examine 
the double layer has led to significant advances in the 
field, mainly due to its ability to explain a large variety 
of colloidal phenomena, and, therefore, has established 
it as the standard representation of the EDL. In turn, 
this adequacy suggests that the RPM already contains 
most of the fundamental traits of a colloidal suspension 
at the usual conditions of experiments and applications. 
However, we consider that the lack of interest in up- 
grading the model of a double layer so as to incorporate 
the effect of ionic size asymmetry stems not only from 
the operational advantages and/or from the past success 
of the RPM but also has been influenced by the com- 
mon belief in the dominance of the counterions in the 
EDL. Such credence has its probable origin in a pioneer- 
ing Poisson-Boltzmann (PB) study by Valleau and Torrie 
24j, where they stated the following: "...we expect the 
double layer properties of a dilute (size-asymmetric) elec- 
trolyte to become similar to those of a completely sym- 
metric electrolyte having an effective size equal to that of 
the counterion. (This remark will be asymptotically ex- 
act for large fields in the Poisson-Boltzmann theory)...". 
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To be more explicit, in Ref. 24] a size-asymmetric elec- 
trolyte next to an electrified wall was analyzed via a quasi 
point-like ions theory known as unequal-radius modified 
Gouy-Chapman (URMGC), which in essence is equiva- 
lent to the classical nonlinear PB equation for a binary 
mixture of punctual ions but with the assignment of dif- 
ferent distances of closest approach (with respect to the 
plate) to anions and cations. Therefore, the remark of 
counterion dominance, quoted above, was indeed formu- 
lated and proved strictly at the Poisson-Boltzmann level. 

Notably, during the past years, a great deal of modern 
treatments of the EDL have also endorsed the idea that 
counterfoils control the properties of hi gh char ged sur- 
faces immersed in both size-symmetric [22], 1371 138( and 
size-asymmetric [13, [H, HH, |32[ primitive model (RPM 
and PM, respectively) electrolytes. Such agreement, 
about the dominance of counterions, between the current 
EDL theories and the old URMGC picture is remarkable 
since in all these new integral equations [13, [H, HU, H3] , 
density functional [H, [32j and mean electrostatic po- 
tential [38j papers the fundamental hypothesis of point- 
ions has been surpassed by including explicitly the hard- 
core and electrostatic correlations neglected in the PB 
theory. Notwithstanding, it must be noted that in the 
cited beyond-PB surveys the "confirmation" of the lead- 
ing role of counterions has been based in studies of 
either charge-asymmetric RPM electrical double layers 
at low/moderate surface charges [13] or, else, of size- 
asymmetric systems near the point of zero charge (PZC) 
[13, HH, H3] ■ In other words, therein the original conclu- 
sion of the preponderance of counterions in the EDL at 
high electric fields has not been really tested. 

In this context, very recently, some of the present 
authors have reported the first theoretical investigation 
of the size-asymmetric spherical electrical double layer 
(SEDL) [lH, where it was found that, contrary to the ac- 
cepted common opinion, for large macrosphere's charge 
densities the counterions do not dominate. As a matter 
of fact, coions are so important that their size can induce 
drastic correlations that bring forth considerable changes 
in the EDL's potential-charge relationship and the surge 
of the charge reversal phenomenon in monovalent salts. 
Remarkably, in the same Ref. [33], the correctness of 
the novel hypernetted chain/mean spherical approxima- 
tion (HNC/MSA) account of the size-asymmetric SEDL 
was already verified, at the level of the radial distribution 
functions (RDFs) , after comparing favorably a few Monte 
Carlo and molecular-dynamics simulations of the ionic 
density profiles with the corresponding HNC/MSA inte- 
gral equation results. Nevertheless, even if this positive 
checking of the HNC/MSA RDFs foresees that other en- 
suing theoretical predictions (e.g. the non-dominance of 
the counterions and the anomalous behavior of the elec- 
trostatic potential at the outer Helmholtz plane) could be 
true, it would be beneficial to have a specific and more 
exhaustive delving of these new features by means of re- 
fined computer "experiments" and/or alternative theo- 
ries (i.e. integral equations, density functionals or mean 



electrostatic potential schemes). Precisely, the primary 
objective of this communication is to extend the research 
of the size-asymmetric SEDL of Ref. [H| by providing 
fresh and comprehensive simulational and theoretical in- 
formation that corroborates the enhancement of the neu- 
tralization and the screening previously found by the the- 
ory and, principally, the non-dominance of counterions at 
high colloidal charges. 

The structure of this paper is as follows: the molecu- 
lar model of the SEDL, theories and Monte Carlo sim- 
ulations are described in Sec. II. In Sec. Ill the partial 
effects of the electrolytic size asymmetry in the nonlinear 
Poisson-Boltzmann scheme via, the URMGC approach 
are discussed, in order to establish a stand point to com- 
pare and discuss the role of the ionic size and size asym- 
metry when these features are included consistently in 
the MC simulations and in the integral equations theory, 
in the HNC/MSA and HNC/HNC approximations. To 
end, a summary of the main findings and some conclud- 
ing remarks are given in Sec. IV. 

II. MODEL AND THEORY 
A. Molecular model 

The main results of this paper are based on the fol- 
lowing representation of the spherical electrical double 
layer (SEDL), which is constituted by a rigid, charged 
spherical colloid of diameter D and surface charge den- 
sity Co , surrounded by a continuum solvent of dielectric 
constant e. The macroion is in contact with two ionic 
species, which are treated as hard spheres of diameters 
Ri (i = 1, 2) with embedded point charges of valences z% 
at their centers. Without loss of generality, we consider 
that i?2 > R\- The interaction potential between the 
macroion, M, and an ion of type i is then given by 

f oo, r<^, 

I e tr ' — 2 > 

where e is the protonic charge. In turn, the interaction 
potential between two ions of species i and j is given by 

Ri+Rj 
oo, r < — 

r> «^ (2) 

iireoer ' — 2 

In the classic literature, the Stern layer or, more prop- 
erly, the Helmholtz surface is the geometrical place cor- 
responding to the closest approach distance between the 
electrolyte ions and the colloid. If we consider an elec- 
trolyte formed by a pair of ionic species of unequal 
size, the inner Helmholtz plane (IHP) is determined by 
(D + Ri)/2, i.e. by the closest approach distance of the 
smallest component to the surface, whereas the outer 
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Helmholtz plane (OHP) is established by (D + R 2 )/2, 
which corresponds to the distance of closest approach for 
the largest species. In the limit of identical sizes the 
IHP and OHP coincide and the usual definition of the 
Helmholtz plane is recovered. 



B. The HNC/HNC and HNC/MSA integral 
equations for the PM-SEDL 

The structural properties of the electrical double layer 
can be obtained from the Ornstein-Zernike equation for 
a multicomponent mixture of S species, which is: 



s 

hij(ri 2 ) = Cij(r u) + y^p; 

1=1 



hii(r 13 )cij(r 32 )dV. (3) 



The set of equations ((3]) requires a second relation, or 
closure, for the functions hij(r) and cy(r). For charged 
systems, the HNC and MSA closures are widely used. 
These relations are given by 



<kj(ri2) = -(3Uij(ru) + h i:j (r 12 ) - ln(h i:j (r 12 ) + 1), (4) 
for HNC, and 



As an especial case of HNC/MSA, if R x = R 2 = in 
Eq. ©, and R 1 ^ and R 2 ^ in Eq. Q, this equation 
reduces to the integral equation version of the nonlinear 
Poisson-Boltzmann (PB) equation (33|: 



I 2 

V 2 ^ = \^ Ztepiexp^-Zteip/kBT). 



(7) 



When Ri and R 2 are equal in Eq. (fTJ), Eq. |(7J) is the 
so-called the Modified Gouy Chapman equation (MGC). 
On the contrary, if Ri and R 2 are different in Eq. (fT|), 
the nonlinear PB equation corresponds to the unequal- 
radius MGC (URMGC) equation 0]. Notice that in 



these quasi point-like ions theories the ionic size is only 
taken into account partially at the level of a closest ap- 
proach distance between the macroion and the ions. Con- 
trastingly, in the MC simulations and in the HNC/MSA 
and HNC /HNC integral equations the ionic size and size 
asymmetry are incorporated consistently, via the primi- 
tive model. 

To better understand and characterize the behavior of 
the SEDL in presence of the electrolytic size asymmetry, 
from the radial distribution functions, gij(r), it is pos- 
sible to calculate two important quantities, namely, the 
integrated charge (IC) 



for MSA, with 



Cij(ri 2 ) = -f3Uij{n 2 ), 



'-./ 1.2....S. 



(5) 



These expressions are complemented by the exact con- 
dition hij = — 1 when < Rij, such as R+j = (Ri + 
^•)/2- 

Let us consider that the species S corresponds to 
macroions at infinite dilution in a binary electrolyte. 
Then Eqs. ([3]) for species S = M and j can be writ- 
ten as: 



2 

hMj(ru) = CMj(n 2 ) + }^pi 
1=1 

J = E2. 



h M i(ri 3 )cij(r 32 )dV. (6) 



When Eq. (fj| is used in © for both CMj(r) and 
Cjj(r) the HNC/HNC integral equation is obtained for 
the SEDL. Besides, if Eq. Q is employed in Eq. © 
only for CMj (f) , arid the cy (r) are approximated by the 
corresponding MSA bulk expressions (i.e. Eq. ([5]) is in- 
serted in Eq. (J3J) for the two electrolytic species), the 
HNC/MSA integral equation is established. The details 
of these i nteg ral equations formalisms can be consulted 
elsewhere [I7ll33l.l3g| and will not be repeated here. How- 
ever, it is important to mention that both schemes satisfy 
the global electroneutrality condition. 



2 pr 

P{r) = z M + V / z lPigi (r)4:iTr 2 dr, (8) 



and the mean electrostatic potential (MEP) 



iP(r) 



t 2 



dt. 



0) 



When the MEP is evaluated at r = (D + R 1 )/2 Eq. © 
corresponds to the MEP at the IHP, which we denote 
as ipiHP- On the other side, if Eq. is calculated at 
r = (D + R 2 )/2 the MEP at the OHP, ipoHP, is obtained. 

With respect to the IC, this quantity is a measure of 
the total net charge inside a sphere of radius r centered 
in the macroion. Then, if D/2 < r < (D + R x )/2 the 
IC is equal to zm, whereas for r — > 00 this quantity 
goes to zero because of the electroneutrality condition. 
Furthermore, the IC has also the property of indicating 
charge reversal when P{t)zm < 0. 



C. Numerical simulations 

Monte Carlo simulations of the SEDL were performed 
considering a cubic box with a macroion fixed at the cen- 
ter under periodic boundaries. Due to the electroneutral- 
ity condition the following relation was satisfied: 



7V_z_ +N+z+ + z M = 0, 



(10) 
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where _/V_ and z_ are the number of ions and the valence 
of the negative species, respectively, N + and z + are the 
number of ions and the valence corresponding to the posi- 
tive species, and zm is the valence of the macroion, which 
is related to the surface charge density as <tq = zmg/t^D 2 . 

In order to take into account the long range nature 
of the coulombic potential, the Ewald sums scheme was 
adopted, using conducting boundary conditions [iol. [4l|. 
The damping constant a was set to a = 5/L and the 
/c-vectors employed to compute the reciprocal space con- 
tribution to the energy satisfied the condition k < 5. The 
length L of the simulation box was assigned considering 
a total number of ions N t = iV_ + N+ « 1000. After 
N t attempts to move an arbitrary ion a Monte Carlo cy- 
cle is counted. The thermalization process consisted of 
2 x 10 4 MC cycles, and from 2 x 10 6 (for high zm values) 
to 6 x 10 6 (for low zm values) MC cycles were completed 
in order to calculate the canonical average. The quality 
of the simulation was tested calculating the IC, which in 
a region far from the macroion and near of the borders 
of the simulation box vanished in all cases, as expected. 





FIG. 1: (Color online:) Radial distribution functions, in- 
tegrated charge and mean electrostatic potential of a size- 
symmetric and size-asymmetric 1:1 salt around a charged 
macroion of valence zm = 4 (ctq = 0.05 C/m 2 ) and diameter 
D — 20 A in the PB approach. The solid and dashed lines 
correspond to URMGC and MGC equations, respectively. 



FIG. 2: (Color online:) The same as in Fig. [Tjbut for zm = 16 
(<to = 0.2 C/m 2 ). 



III. RESULTS AND DISCUSSION 

A. The role of the ionic size asymmetry in the PB 
scheme 

Physically, in the MGC and URMGC equations the 
ionic size is considered only partially because the elec- 
trolytic ions are allowed to be close to the macroion un- 
til a closest approach distance for each species, but the 
ions interact among them as charged points. Thus, the 
MGC and URMGC equations correspond, to the simplest 
manner in which the effects of the ionic size and the size 
asymmetry can be studied in the EDL, as it was done 
by Valleau et al. in the 1980s for the planar instance 
24]. Notwithstanding, although URMGC represented a 
step beyond MGC, when the ionic size and in particu- 
lar the ionic size asymmetry is fully taken into account 
(as occurs in the MC simulations and in the HNC/MSA 
and HNC/HNC integral equations) new features absent 
in PB picture emerge. Therefore, in order to later com- 
pare and discuss the consequences of a complete con- 
sideration of ionic size asymmetry in the SEDL, we will 
begin with a review of size asymmetry in the URMGC 
approach, where the excluded volume effects are embod- 
ied partially only in the colloid-ion interactions. Thus, 
let us first study a spherical macroion of diameter and 
surface charge density D = 20 A and er = zmg/ttD 2 , 
respectively , surrounded by a 1:1, 1 M electrolyte, in 
a continuum solvent of dielectric constant e = 78.5 at 
a temperature T = 298 K . In the size-symmetric case 
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FIG. 3: (Color online:) Mean electrostatic potential at the 
IHP as function of the surface charge density cro, for a 
1:1, 1 M electrolyte around a macroion of diameter D — 
20 A. The simulations results were calculated for zm = 
0, 4, 8, 12, 16, 20, 24, 28 and 32. The diameters of the ionic 
species in the PM are JJ_ = 4.25 A and R+ = 8.5 A. In 
the RPM, the ionic diameters are equal to the counterions 
in the PM, i.e. 4.25 A . The maximum approach distances 
dij for ion-ion and macroion-ion interactions for theory and 
simulation are given by Eqs. (11111 y (|12|l . The triangles and 
the circles correspond to Monte Carlo simulation results in 
the PM, MCpa/, and in the RPM, MCrpm, respectively. 
The dotted and solid lines correspond to HNC/MSApm and 
HNC/MSA_rpm, and the dot-dashed and dashed lines are as- 
sociated to HNC/HNCpm and HNC/HNChpm, respectively. 
The dashed line with multiplication symbols denotes URMGC 
and the dotted line with plus symbols is for MGC. 



(i.e. for the MGC theory), the maximum approach dis- 
tance for both species is 4.25/2 A , whereas in the size- 
asymmetric (i.e. for the URMGC theory) is 4.25/2 A for 
anions and 8.5/2 A for cations. Since we will consider 
only cro > values, in both instances the counterions 
have the same properties, being the size of the coions the 
unique difference between the MGC and URMGC sys- 
tems. Therefore, in Fig. [T] we compare the RDFs, ICs 
and MEPs curves associated to the size-symmetric and 
to the size-asymmetric cases, when the valence of the 
macroion is zm = 4 (cro = 0.05C/m 2 ). At the level of 
the RDFs, in Fig. QJ, it is observed that the MGC profiles 
for the size-symmetric case enclose the RDFs of the size- 
asymmetric electrolyte described by URMGC. Besides, 
from Figs. [TJd andQJ: it is seen that the region not allowed 
for big cations but accessible for the small anions in the 



0.1 0.2 0.3 0.4 

% (C/m 2 ) 

FIG. 4: (Color online:) The same as in Fig. [3] but for a 2:2 
electrolyte at 0.5 M concentration. 



URMGC theory contributes significantly to the increase 
of the neutralization and the screening of the spherical 
EDL when contrasted with the MGC theory (notice that 

P(r)uRMGC < P{t)mgc and ip{r) URM GC < 4>{r)MGC 
for all the r plotted). These differences in the RDFs, 
P(r), and ip{ r ) are expected to augment if zm decreases, 
with the largest dissimilarities occurring precisely at the 
point of zero charge (PZC). Additionally, it is foreseen 
a MEP equal to zero at the closest approach distance 
of anions in the MGC results and an electrostatic poten- 
tial different from zero for URMGC at the same distance. 
This happens for a z : z salt when zm = because, under 
these conditions, the RDFs for MGC must coincide for 
symmetry reasons, while for URMGC the ionic density 
profiles must be completely different in order to fulfill the 
electroneutrality condition (i.e. as the macroion is un- 
charged, only the accumulation of the large ionic species 
after the OHP can compensate the adsorbed charge of 
the small ions in between the Helmholtz planes). On the 
contrary, if zm increases and there are no crossings be- 
tween the RDFs of counterions and coions for MGC, and 
the same behavior is displayed by URMGC, the MGC 
curves are the limit of URMGC profiles when zm —> oo 
due again to the electroneutrality condition. Such phe- 
nomenon has been already discussed in (33[. Thus, even 
if, strictly, the MGC and URMGC profiles should not be 
the same for high zm values (because there is always a 
region where the small coions can exist in URMGC, see 
FiglSk); the structural properties of the EDL, as the in- 
tegrated charge and the mean electrostatic potential, are 
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FIG. 5: (Color online:) Mean electrostatic potential at the 
OHP as function of the surface charge density <tq, around a 
macroion of diameter D = 20 A. The simulation results were 
calculated for z M = 0,4,8,12,16,20,24,28 and 32. The di- 
ameters of the ionic species in the PM are R- — 4.25 A and 
R+ = 8.5 A. The maximum approach distances dij for ion-ion 
and macroion-ion interactions for theory and simulation are 
given by Eqs. (|lip y (112[) . The triangles correspond to Monte 
Carlo simulations. The dotted, dot-dashed and dashed are as- 
sociated to HNC/MSApm and HNC/HNCpm, and URMGC, 
respectively. In Fig. [5^ the electrolyte is 1:1, 1 M, whereas in 
Fig. [5)3 the salt is 2:2, 0.5 M. 
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FIG. 6: (Color online:) Integrated charge and mean elec- 
trostatic potential as function of the distance for a 1:1, 1 
M electrolyte near a macroion of diameter D = 20 A and 
2jvi = 24. The diameters of the ionic species in the PM are 
R- = 4.25 A and R+ = 8.5 A. In the RPM, the ionic diame- 
ters are equal to the counterions in the PM, i.e. 4.25 A. The 
maximum approach distances dij for ion-ion and macroion-ion 
interactions for theory and simulation are given by Eqs. (|11[1 
y (|12[) . The triangles and the circles correspond to MCpm 
and MC rpm, respectively. The dotted and solid lines cor- 
respond to HNC/MSApA f y HNC/MSAppm, and the dot- 
dashed and dashed lines are associated to HNC/HNCpm and 
HNC/HNCppa/, respectively. The dashed line with multi- 
plication symbols denotes URMGC and the dotted line with 
plus symbols is for MGC. The distance r' is measured from 
the macroion's surface. 



indeed very similar as can be observed in Figs. [2]d and 
[2J;. Consequently, far from the PZC the properties of the 
spherical EDL are expected to be practically the same 
for URMGC and MGC because the coion's contribution 
is negligible and the counterions are the same in both 
cases. This last merging between URMGC and MGC is 
precisely the so-called dominance of counterions in the 
spherical EDL and will be of decisive importance in the 



corresponding potential-charge relationship, as it will be 
shown later. 
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FIG. 7: (Color online:) The same as in Fig. [6] but for a 2:2, 
0.5 M salt. 



B. The role of the ionic size asymmetry in the 
primitive model: MC simulations and theoretical 
results 



As discussed, one of the main differences between size- 
symmetric and size-asymmetric EDLs salt in the PB 
viewpoint is the increment in the neutralization or screen- 
ing predicted by URMGC, with respect to MGC, near 
the PZC. In addition, for high surface charge densities 
the properties of the SEDL are expected to be practically 
the same, such as the dominance of counterions arises in 
a quasi-point like description (i.e. MGC and URMGC). 

Now, in this section we will show comprehensive MC 
data and integral equations (IE) results in which the ionic 
size asymmetry is taken into account consistently in the 
Eqs. fl]) and @ (and not only at the macroion-ion level 
as it was done in the nonlinear PB equation of Sec. Ill A) 
in order to display the effects in the SEDLs due to a more 
realistic treatment of size-symmetric and size-asymmetric 
salts. 



In all the following simulations and theoretical calcu- 
lations we considered a macroion of diameter D = 20 
A and ao > 0, immersed in a continuum solvent of di- 
electric constant e = 78.5 at a temperature T = 298 
K, in presence of a binary electrolyte. In the primitive 
model (PM) the diameter of the counterions is R- — 4.25 
A and the diameter of coions is the double, i.e. R + = 8.5 
A. For the restricted primitive model, the size of both 
species is the same of the counterions in the PM, i.e. 
4.25 A. Thus, the maximum approach distances in the 
PM and the RPM correspond to those of the electrolyte 
for URMGC and MCG, respectively. This information 
is summarized in Eqs. (fTTj) and (fT2|) . Notice that the 
diameter of the macroion and the ionic size asymmetry 
correspond to values that emphasize the spherical geom- 
etry of the EDL and that have been typically used in 
previous works [U H S El ■ 



Very importantly, the mean electrostatic potential at 
some distance from the surface of the macroion is fre- 
quently associated with the so-called electrokinetic 
tential at the shear plane (or the zeta potential, £) 
Such quantity is very relevant in colloidal studies since it 
is experimentally measurable and allows to characterize 
and summarize the behavior of the SEDL, as functions 
of the colloidal charge, in a single potential-charge plot. 
Furthermore, the zeta potential is used often in physical- 
chemistry to characterize the macroscopic properties and 
the stability of charged colloidal dispersions [5j, l(| . Given 
that, in the past, the ipmp or the ipoHP have identified 
with £, we start by showing the ip-o~o curves at the IHP 
and at the OHP for our PM systems. In Fig. H MC sim- 
ulations of the mean electrostatic potential at the IHP for 
a 1 M, 1:1 electrolyte in the PM and the RPM are shown. 
The first notable feature there is the merging of the MGC 
and URMGC curves for high ao values. Precisely, this 
asymptotic conduct illustrates the dominance of counte- 
rions at the level of tp-ao- Besides, in this figure it is seen 
that the maximum difference between the MEPs corre- 
sponding to MGC and URMGC happens precisely at the 
PZC, as it had been pointed out in Sec. Ill A. In strong 
contrast, the most evident characteristic displayed by the 
simulational data of ifjiHP for the size-symmetric (RPM) 
and size-asymmetric (PM) instances is that these curves 
do not converge to the same one when ao augments. This 
confirms that the counterions do not always dominate the 
EDL far from the PZC, or, in other words, exemplify 
the importance of the size of the coions at high colloidal 
charges, as it was theoretically predicted in Ref. j3$i] . Ad- 
ditionally, the simulational data shows a potential differ- 
ent from zero at the PZC for the size-asymmetric elec- 
trolyte. Clearly, such behavior is due to the fact that the 
small negative ions are allowed to be closer to the surface 
than the big positive ions, i.e. for ao — the negative sign 
of the MEP at the IHP results from the size asymmetry 
of the 1:1 electrolyte. Interestingly, analogous results had 
been theoretically predicted for the RPM planar EDL of 
charge asymmetric species [3^,13]. However, these data 
have not been confirmed simulationally. With respect 
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d- = and d+ = for MCpm, HNC/MSA pM 

HNC/HNC PM and URMGC; 

d _= d+ = £+£^ f or MCrpm, HNC/MSA ppM 

HNC/HNC PPM and MGC; 




d_ _ = R. and d++=R+, for MCpm , HNC/MSA PM , HNC/HNC PM ; 

d-+ = d+- = for MCpm , HNC/MSA PM , HNC/HNC PM ; 

d_ _ = d ++ = d_+ = d+- = for MCppm , HNC/MSA PPM and HNC/HNC 

d- - = d+ + = d- + = d+ - = 0, for URMGC and MGC. 



to the performance of the integral equation theories in 
the RPM, it is remarkable that both the HNC/HNC and 
HNC/MSA theories agree with the simulation data, al- 
though, for the PM case, HNC/HNC follows closely the 
MC curves than HNC/MSA. 

In addition, in Fig. [3] the MC simulations and in- 
tegral equations (IE) predict a larger value of <to for 
which the negative sign of the MEP at the IHP changes 
to positive than that corresponding to URMGC. This 
exemplify again the importance of the entropic contri- 
butions in the adsorption of counterions between the 
Helmholtz planes when the ionic size correlations are 
considered consistently (as occurs in the MC simula- 
tions), in contrast with its partial inclusion, when only 
different approach distances in the macroion-ion interac- 
tion are considered, which neglects ionic excluded vol- 
ume effects outside of the OHP (as in the URMGC the- 
ory). Another phenomenon observed in Fig. [3] is that 
-0(cto)pm < '0(c r o)i?PAf at the IHP for all <7q values plot- 
ted. Near the point of zero charge this is explained, in 
the PM, in terms of the adsorption of negative coun- 
terions that are not neutralized by the big coions, as 
it happens in the RPM case. When (To augments, the 
positive bare charge overcomes the contribution of the 
negative counterions near the macroion's surface and the 
MEP's sign changes from negative to positive as it was 
mentioned above. Nevertheless, for all the <tq values 
displayed the counterions in the PM provide an extra 
screening which is not present in the RPM, which leads 
to the ip(<Jo)pM < "0( (T o)fiPA/ condition. Moreover, in 
the monovalent case of the PM this extra screening can 
be related not only with a higher adsorption of negative 
counterions with respect to the RPM but also with the 
presence of charge reversal, that is absent in the RPM. 
This behavior will be clarified later when the correspond- 
ing P(r) profiles be presented. 

With the purpose of performing a more stringent test 
for the theories, in Fig. [4] we present MC simulations 
of the MEP at the IHP for a 0.5 M, 2:2 electrolyte in 
the PM and the RPM. In this instance, the features al- 
ready observed in the simulations of monovalent ions are 
accentuated. In particular, the most important finding 
is the corroboration of the non-dominance of counteri- 
ons for divalent ions. In addition, in the size-asymmetric 
case a very strong adsorption of negative counterions in 



the PM is also observed. The importance of excluded 
volume effects is evinced by noticing that near the PZC 
the interval of oo for which the ipmp is negative is larger 
for MC simulations and IE theories than for URMGC. 
Furthermore, when <7o increases after ip(< 7 o)pM > 0, the 
ipiHP reaches a maximum and for still larger do values 
ip(ao)pM < again. The appearance of a maximum 
in the ip( a o) pl°t f° r the RPM is related to presence of 
charge reversal. Besides, the early MEP's change of sign 
in the PM after the maximum displayed by MC simu- 
lations suggests an extra adsorption of counterions with 
respect to the RPM case, i.e. it is expected an accentu- 
ated charge reversal that screens more strongly the posi- 
tive bare charge of the macroion for high <jq values. This 
behavior will be clearly exhibited in the corresponding 
P(r) profiles later. On the other hand, the simulational 
confirmation of the reentrance in the sign of the MEP for 
divalent ions in the PM, i.e. the double change of sign 
of ipiHP, suggests the possibility of observing a corre- 
sponding reentrance in the experimental electrophoretic 
mobility (/x), if the Smolouchowski equation (/i = eC/v) is 
valid, as it had been theoretically foreseen by HNC/MSA 
[33| for a larger macroion. This means that the ionic size 
asymmetry could then cause a reversed mobility in the 
motion of a macroion in an electrophoresis experiment 
near the PZC, changing to the "correct" direction when 
its surface charge density augments, but inverting again 
its movement at high 00 values. Also, notice that for 
the size-asymmetric instance HNC/HNC shows a better 
agreement with the simulation data than HNC/MSA fol- 
low colloidal charges (when <tq < 0.16 C/m 2 approxi- 
mately). Contrastingly, for high surface charge densities 
((To < 0.16 C/m 2 ) the opposite behavior is observed, i.e. 
a substantial deviation from the MC data is displayed by 
HNC/HNC, which contrasts with the good accordance 
shown by HNC/MSA. This conduct is remarkable since, 
for high colloidal charges and high entropic-electrostatic 
ionic correlations, in the divalent case HNC/MSA is bet- 
ter than HNC/HNC in opposition to the univalent in- 
stance. 

In Fig. [5] the simulational results of the MEP at the 
OHP, ipoHP, for the PM monovalent and divalent salts 
are plotted. As we said before, the importance of the 
study of this curves arises from the fact that the ipoHP 
could be identified with the zeta potential (£). In Fig. 
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[5^,, which corresponds to the 1:1, 1 M electrolyte, a non- 
monotonic behavior of the MEP as a function of a is 
observed. This conduct is reproduced correctly by IE the- 
ories, with HNC/HNC showing a better agreement than 
HNC/MSA. Contrastingly, URMGC presents a mono- 
tonic behavior, which completely differs from the simu- 
lation data, and predicts large values of the electrostatic 
potential for high o~o as it can be seen in the inset. The 
simulational results of ipoHP for the 2:2, 0.5 M electrolyte 
are portrayed in Fig. [SJd. Here, in the MC simulations 
it is observed that for any value of cto the ipoHP is nega- 
tive and decreases monotonically as a function of a®. A 
similar curvature has been theoretically predicted in the 
RPM spherical EDL by one of the present authors for 
a macroion of diameter D = 80 A immersed in a 2:2, 
0.5 M electrolyte of ionic diameter equal to 7 A [2l|, at 
approximately the same volume fraction used here in our 
PM. Thus, the simulational results presented in this work 
corroborate that such anomalous curvatures in the MEP, 
at high ionic volume fractions, are a real feature in the 
primitive model. Even if this phenomenon is interesting 
just from the theoretical point of view [H, could 
also be relevant in the description of non-intuitive at- 
tributes of double layer systems such as the occurrence of 
negative differential capacitances [Z^, 0, IH, [EH, [13] . 
Once again, the whole behavior displayed by the simu- 
lations is well captured by the the IE theories, although 
now HNC/MSA is closer to simulations than HNC/HNC. 
Contrastingly, URMGC exhibits a monotonic behavior in 
which the MEP increases as a function of ao, as occurred 
in the 1:1 instance. 

As was noticed for the ip(°~o) relationship, one of the 
consequences of the ionic size asymmetry in the PM elec- 
trical double layer, is the enhancement of the neutraliza- 
tion and screening at high surface charges. To illustrate 
this in terms of the ionic charge adsorption, in Fig. [5^ we 
have plotted the integrated charge profiles for the 1:1, 1 
M electrolyte in the PM and RPM, when the macroion's 
valence is Zm — 24. Here, it is clearly seen that the ionic 
size asymmetry not only promotes a higher adsorption of 
counterions, i.e. P(r)pM < P(v)rpm, but also that can 
induce the appearance of charge reversal in monovalent 
electrolytes, showing a minimum at r'/i?_ « 2. Con- 
sequently, MC simulations hint that charge reversal can 
occur even in presence of monovalent salts whenever the 
high coupling conditions are present, i.e. high electrolyte 
concentration or large hydration of the electrolyte. Note 
that the IE theories reproduce very well the charge rever- 
sal behavior, whereas an incorrect monotonic decrease of 
the IC is shown by URMGC. 

The overcompensation of the native charge in the PM 
is also reflected in the corresponding MEP curves as a 
function of the distance as can be verified in Fig. [5Jd. 
In this case the monotonicity of the ip(r)ppM and the 
non-monotonic behavior of ip{r) pm can be easily deduced 
from the corresponding P{r) profiles and the Eq. ©. In 
particular, note that the condition P(r)pM < P{v)rpm 
implies that tp(r)pM < V'( r )-RPM for MC simulations and 



HNC/HNC and HNC/MSA theories, i.e. to a higher neu- 
tralization of the macroion's bare charge (as observed in 
P(r)) a higher screening in ip( r ) is associated. 

In the Fig. [7^ the charge adsorption for a 2:2, 0.5 
M electrolyte as function of the distance in the PM and 
RPM is displayed when the valence of the macroion is 
zm = 24. Notice that the integrated charge simulation 
curve shows that in the RPM electrical double layer there 
is charge reversal, contrasting with the behavior of the 
1 M monovalent case previously portrayed, where it was 
found the absence of this feature even for a higher elec- 
trolytic concentration. Besides, when the ionic size asym- 
metry is present in the PM the charge reversal, already 
observed in the RPM, is notably enhanced. This illus- 
trates the fact that the high electrostatic- entropic cou- 
pling conditions for which the charge reversal appears can 
be relaxed for multivalent salts, i.e., it is expected that for 
1 : z salts the ionic size and/or the ionic size asymme- 
try (coming from the ionic hydration for example) be- 
come very important even at moderate salt concentra- 
tions reachable experimentally, as it has been rep orted 
by several experimental works [H, H3, [H, [56|, HJ. The 
MEP curves corresponding to the IC profiles discussed 
previously now are plotted in Fig. [TJd. Consistently, the 
MC simulations and the HNC/HNC HNC/MSA theories 
predict that near the macroion's surface the screening 
in the PM is higher than in the RPM, with the MEP 
presenting a non-monotonic behavior in both cases. Fur- 
thermore, here it is noticed that the overestimation of 
the screening in the HNC/HNC potential profile has its 
origin in the charge overcompensation displayed by the 
corresponding P(r) graphed in Fig. [7Ji. 



IV. CONCLUDING REMARKS 

Monte Carlo simulations of the primitive model spher- 
ical electrical double layer, in the presence of cither 
monovalent or divalent salts, were performed and com- 
pared with data of the HNC/HNC and HNC/MSA inte- 
gral equations and of the the MGC and URMGC quasi- 
punctual-ions schemes. One of the most simple man- 
ners in which the ionic size and size asymmetry can be 
taken into account partially in the EDL is via the MGC 
and URMGC approaches, in which the excluded volume 
effects are considered solely in the macroion-ion inter- 
actions by means of a closest approach distance. When 
size-symmetric and size-asymmetric semi-punctual EDLs 
(i.e. systems with either equal or different closest ap- 
proach distances for the otherwise-punctual electrolytic 
species) with the same type of counterions are consid- 
ered by the MGC and URMGC formalisms, respectively, 
it was exhibited here that, at low colloidal charges, one of 
the main effects of including the ionic size asymmetry in 
the URMGC theory is an enhancement of the screening 
and the neutralization in the EDL with respect to the 
MGC (or size-symmetric) instance. On the other hand, 
far from the point of zero charge, it was shown that the 
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RDFs, ICs, and MEPs predicted by the MGC and UR- 
MGC equations displayed always a monotonic comport- 
ment, and that the ionic distributions almost coincided 
for all the values of r, with exception, maybe, of the zone 
comprised between the Hclmholtz planes. However, as 
(7 — ► oo, the diferences between the MGC and URMGC 
radial distribution functions, and therefore between all 
their concomitant structural and thermodynamic prop- 
erties, go asymptotically to zero. This behavior is the so- 
called dominance of the counterions. Contrastingly, when 
the ionic finite size and size asymmetry are embodied 
consistenly into the macroparticle-ion and ion-ion inter- 
actions, as occurred in the Monte Carlo (MC) simulations 
and the HNC/HNC and HNC/MSA integral equations, 
several characteristics absent in the nonlinear Poisson- 
Boltzmann data (e.g. the non-montonic behavior of the 
RDFs, ICs, and MEPs and the phenomenon of charge re- 
versal) emerge. Most importantly, our MC simulations of 
the primitive model EDL corroborated the fact that the 
ionic size asymmetry augments the colloidal charge neu- 
tralization and the screening in comparison with the size- 
symmetric case even at high values of 00 , proving the non- 
dominance of the counterions in the primitive model. On 
the theoretical side, the HNC/HNC and HNC/MSA inte- 
gral equations displayed consistent results with MC simu- 
lations, without a notable predominance in the accuracy 
between them, whereas MGC and URMGC evidenced 
the limitations of the Poisson-Boltzmann theories. Other 
consequences of the ionic size asymmetry in the EDL 
that have been confirmed by the present numerical simu- 
lations and HNC/HNC and HNC/MSA calculations are 
the appearance of charge reversal in monovalent salts and 
the reentrance of the mean electrostatic potential at the 
outer Hclmholtz plane (i.e., the change of sign of the 
4>Ohp from negative to positive and, then, to negative 
again when uq increases from zero) for divalent salts. If 
the usual identification between the well-known electroki- 
netic zeta potential and the mean electrostatic potential 



in the neighborhood of the Helmholtz region is assumed 
0, such reentrance in tpoHP could be of relevance for 
mobility experiments since it indicates the possibility of 
observing an alternating direction in the electrophoresis 
of a colloid, immersed in a multivalent electrolytic bath, 
as a function of the native macroparticle's charge. Fur- 
thermore, the reported MC simulations have also evinced 
that anomalous curvatures can appear at the OHP in 
the primitive model EDL, which could be important for 
several recent investigations about non-intuitive phenom- 
ena (e.gjjthe appearance of negative differential capac- 
itances IH, HH, HI) in electrolyte-electrode 
systems. In summary, the data reported in this paper 
suggest that the ionic size and, especially, the ionic size 
asymmetry should be considered as very sensitive param- 
eters that, in combination with the concentration and va- 
lence of the electrolyte and the macroion's surface charge 
density, control the electrostatic-entropic coupling in the 
primitive model EDL. In particular, given that for mul- 
tivalent salts the ionic hydration augments notably the 
finite size and size asymmetry effects, this could represent 
a way to attain a high electrostatic-entropic coupling at 
reasonable experimental conditions for which the charge 
reversal could be detected [HI, [H, [H, HE Hz| • 
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